home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / ARGONET / PD / MATHS / RLAB / RLAB125.ZIP / !RLaB / toolbox / erf < prev    next >
Text File  |  1995-11-20  |  2KB  |  107 lines
  1. //-------------------------------------------------------------------//
  2. //
  3. // Syntax:    erf    ( A )
  4. //        erfc   ( A )
  5. //        erfcc  ( A )
  6. //        inverf ( A )
  7.  
  8. // Description
  9.  
  10. // Error function, complimentary error functions, and inverse error
  11. // function. 
  12.  
  13. // erf   (A) = 2/sqrt(pi) sum exp(-t^2) (from 0 to A)
  14. // erfc  (A) = 1 - erf (A)
  15. // erfcc (A) = erfc (A) (but faster)
  16.  
  17. // See Also: gamma
  18.  
  19. //-------------------------------------------------------------------//
  20.  
  21. // Dependencies
  22.  
  23. require gamma
  24.  
  25. erf = function ( X )
  26. {
  27.   local (e, i);
  28.  
  29.   e = zeros (size (X));
  30.   for (i in 1:X.n)
  31.   {
  32.     if (X[i] < 0)
  33.     {
  34.       e[i] = -gammp (X[i].^2, 0.5);
  35.     else
  36.       e[i] =  gammp (X[i].^2, 0.5);
  37.     }
  38.   }
  39.   return e;
  40. };
  41.  
  42. erfc = function ( X )
  43. {
  44.   local (e, i);
  45.  
  46.   e = zeros (size (X));
  47.   for (i in 1:X.n)
  48.   {
  49.     if (X[i] < 0)
  50.     {
  51.       e[i] = 1.0 + gammp (X[i].^2, 0.5);
  52.     else
  53.       e[i] = gammq (X[i].^2, 0.5);
  54.     }
  55.   }
  56.   return e;
  57. };
  58.  
  59. //
  60. // A simpler, and much faster version of erfc().
  61. //
  62.  
  63. erfcc = function ( X )
  64. {
  65.   local (ans, i, t, z);
  66.  
  67.   z = abs (X);
  68.   t = 1 ./ (1 + 0.5*z);
  69.   ans = t.*exp (-z.*z-1.26551223+t.*(+1.00002368+t.*(0.37409196+...
  70.                  t.*(+0.09678418+t.*(-0.18628806+t.*(0.27886807+...
  71.                  t.*(-1.13520398+t.*(+1.48851587+...
  72.                  t.*(-0.82215223+t.*0.17087277)))))))));
  73.   for (i in 1:X.n)
  74.   {
  75.     if (X[i] < 0)
  76.     {
  77.       ans[i] =  2-ans[i];
  78.     }
  79.   }
  80.   return ans;
  81. };
  82.  
  83. inverf = function ( X )
  84. {
  85.   local (b, bp, f, fd, m, n, tol, tp);
  86.   global (pi)
  87.  
  88.   //
  89.   // Newton-Raphson iteration on ERF
  90.   //
  91.  
  92.   tol = 0.00001;
  93.   tp = 1/sqrt (2*pi);
  94.   m = X.nr;
  95.   n = X.nc;
  96.   bp = ones (m, n);
  97.   b = zeros (m, n);
  98.   while (any (any (abs (b-bp) > tol)))
  99.   {
  100.     bp = b;
  101.     f = ((1-erfcc (bp)) - X)/2;
  102.     fd = exp (-0.5*bp.*bp).*tp;
  103.     b = bp - f./fd;
  104.   }
  105.   return b;
  106. };
  107.